获取基因有效长度的N种方法

最近有粉丝自告奋勇希望可以把他自己在简书等平台的生物信息学笔记分享在我们生信技能树公众号，在专业的舞台上跟大家切磋！

非常欢迎，他前面的分享是：Counts FPKM RPKM TPM CPM 的转化 

下面是简书的昵称:  嘿嘿嘿嘿哈  的分享

在RNAseq的下游分析中，一般都会将上游处理完得到的原始counts数转变为FPKM/RPKM或是TPM来进行后续的展示或分析，其定义和计算公式在 前面的分享是：Counts FPKM RPKM TPM CPM 的转化 提到了。

需要注意一点的是，在计算FPKM/RPKM和TPM时，基因长度一般指的都是基因有效长度effective length，即该基因的外显子总长度或转录本总长度，以此为标准来消除测序造成的基因长度影响才更为准确。参见生信技能树文章：

基因长度之多少 | 生信菜鸟团 (bio-info-trainee.com)

那么问题来了，在计算FPKM/RPKM时，每个基因的基因有效长度数据该如何获取呢？
我总结了几种获取基因有效长度（或非冗余总外显子长度、总转录本长度）的方法，现整理如下：

一、从上游输出文件结果中获取基因有效长度

一般而言，RNA-seq得到原始counts表达矩阵最常用到的上游软件就是featureCounts和Salmon了，在这两类软件的输出结果中，除了基因（或转录本）的counts信息外，也包含了基因有效长度信息，如featureCounts输出文件的Length这一列对应的就是基因有效长度。（P.S. 之前一直以为featureCounts的Length只是单纯的基因长度，后来经过多种方法比较后发现其实Length这一列就已经是基因的有效长度了...在文章后面我也会展示这几种方法比较的结果）

因此，最方便的做法就是在下游获取counts矩阵时，将基因有效长度信息也同时提取出来用于后续的基因表达量转化。

1. 针对featureCounts的输出文件

在R中读取featureCounts的输出文件，提取Length和对应的geneid信息，再按照counts中的rowname（geneid）匹配排序，即可进行后续的TPM、FPKM值的计算了。

featurecounts输出文件格式

rm(list=ls()) 
options(stringsAsFactors = F)  
library(tidyverse) 
# ggplot2 stringer dplyr tidyr readr purrr  tibble forcats 
library(data.table) #可多核读取文件 
a1 <- fread('counts.txt', header = T, data.table = F)#载入counts，第一列设置为列名 
### counts矩阵的构建 
counts <- a1[,7:ncol(a1)] 
#截取样本基因表达量的counts部分作为counts  
rownames(counts) <- a1$Geneid #将基因名作为行名 
### 从featurecounts 原始输出文件counts.txt中提取Geneid、Length(转录本长度)， 
geneid_efflen <- subset(a1,select = c("Geneid","Length"))        
colnames(geneid_efflen) <- c("geneid","efflen")   
geneid_efflen_fc <- geneid_efflen #用于之后比较 
### 取出counts中geneid的对应的efflen 
dim(geneid_efflen) 
efflen <- geneid_efflen[match(rownames(counts),                               
                              geneid_efflen$geneid),"efflen"] ### 计算 TPM #TPM (Transcripts Per Kilobase Million)  每千个碱基的转录每百万映射读取的Transcripts 
counts2TPM <- function(count=count, efflength=efflen){   RPK <- count/(efflength/1000)  #每千碱基reads (“per million” scaling factor) 长度标准化   
                                                      PMSC_rpk <- sum(RPK)/1e6   #RPK的每百万缩放因子 (“per million” scaling factor ) 深度标准化  
                                                      RPK/PMSC_rpk                      
                                                     }   
tpm <- as.data.frame(apply(counts,2,counts2TPM)) colSums(tpm)

其中geneid_efflen内容如下

geneid_efflen

2. 针对Salmon的输出文件

Salmon的输出结果以及各内容的解释如下。Salmon Output File Formats - Salmon 1.8.0 documentation
值得注意的是，salmon不仅给出了基因有效长度Length（转录本长度），还给出了EffectiveLength，即经过考虑各种因素矫正后的转录本有效长度。官方更推荐使用EffectiveLength进行后续的分析，它结果中的TPM值也是根据EffectiveLength计算的。

Salmon的输出结果

Salmon的输出结果官方解释

我们一般使用tximport导入salmon的输出文件“quant.sf”（转录本的统计结果）和转录本id与gene symbol对应关系文件，会生成以下几个数据：
"abundance" "counts" "length" "countsFromAbundance"
tximport生成的Length就是EffectiveLength，而"abundance" 就是TPM值，我们提取Length用于后续计算FPKM。注意，由于每个样本中基因的EffectiveLength有差异，我们要提取的实际为EffectiveLength的矩阵(或者可以每行EffectiveLength取均值)。

library(tximport)  #t2s为从gtf文件中提取的transcript_id和symbol的对应关系文件 
t2s <- fread("t2s_vm29_gencode.txt", data.table = F, header = F) 

##创建quant.sf所在路径  导入salmon文件处理汇总 
quantdir <- file.path(getwd(),'salmon'); 
quantdir files <- list.files(pattern="*quant.sf",quantdir,recursive=T); files  #显示目录下所有符合要求的文件
files <- file.path(quantdir,files);files 
txi_gene <- tximport(files, type = "salmon", tx2gene = t2s) 

##提取出counts/tpm表达矩阵 
counts <- apply(txi_gene$counts,2,as.integer) #将counts数取整 
rownames(counts) <- rownames(txi_gene$counts)
tpm <- txi_gene$abundance  ###abundance为基因的Tpm值 

###提取geneid_efflen_mat 
geneid_efflen_mat <- txi_gene$length  ###Length为基因的转录本有效长度 

## 计算 TPM 、FPKM  
if (F) { #可直接从txi的"abundance"  中提取，不用运行   
  tpm <- data.frame(rownames(counts),row.names = rownames(counts))     
  for (i in 1:ncol(counts)) {       
    count <- counts[,i]        
    efflength <- geneid_efflen_mat[,i]       
    RPK <- count/(efflength/1000)   #每千碱基reads (reads per million) 长度标准化       
    PMSC_rpk <- sum(RPK)/1e6        #RPK每百万缩放因子 (“per million” scaling factor ) 深度标准化 
    tpm00 <- RPK/PMSC_rpk         
    tpm <- data.frame(tpm,tpm00)      
    rm(tpm00)    
  }    
  tpm <- tpm[,-1];  
  colnames(tpm) <- colnames(counts);  
  head(tpm)       
}      
## 计算 fpkm   
if(T){     
  fpkm <- data.frame(rownames(counts),row.names = rownames(counts))    
  for (i in 1:ncol(counts)) {    
    count <- counts[,i]    
    efflength <- geneid_efflen_mat[,i]   
    PMSC_counts <- sum(count)/1e6   #counts的每百万缩放因子 (“per million” scaling factor) 深度标准化    
    FPM <- count/PMSC_counts        #每百万reads/Fragments (Reads/Fragments Per Million) 长度标准化      
    fpkm00 <- FPM/(efflength/1000)     
    fpkm <- data.frame(fpkm,fpkm00)    
    rm(fpkm00)   
  }     
  fpkm <- fpkm[,-1];
  colnames(fpkm) <- colnames(counts) 
} 

如果想要提取一般意义上的基因有效长度，需要利用“quant.genes.sf”文件（基因的统计结果，需要在进行salmon时加上参数 -g ，后接gtf文件），提取Length这一列的信息。

a2 <- fread("quant.genes.sf",             
            data.table = F) 
geneid_efflen <- subset(a2, select = c("Name","Length"))
colnames(geneid_efflen) <- c("geneid","efflen")  
geneid_efflen_salmon <- geneid_efflen #用于之后比较

二、 从gtf文件中计算获取基因有效长度

整理了两种从gtf文件中计算获取基因有效长度的方法（非冗余外显子长度之和），参考这两篇文章：
基因长度并不是end-start - 简书 (jianshu.com)Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)
由于处理数据量很大，代码速度运行较慢，因此在以下代码中还调用了parallel包进行多核运算处理

1. 利用GenomicFeatures包导入gtf处理

library(parallel) #并行计算  parApply parLapply parSaplly 
cl <- makeCluster(0.75*detectCores())  #设计启用计算机3/4的核

## 利用GenomicFeatures包导入gtf处理
    library(GenomicFeatures)
    txdb <- makeTxDbFromGFF("gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz",
 format="gtf") 
    exons_gene <- exonsBy(txdb, by = "gene") ###提取基因外显子
    head(exons_gene)

    ##计算总外显子长度：用reduce去除掉重叠冗余的部分，,width统计长度，最后计算总长度
    exons_gene_lens <- parLapply(cl,exons_gene,function(x){sum(width(reduce(x)))}) 
    exons_gene_lens[1:10]
    
    ##转换为dataframe
    geneid_efflen <- data.frame(geneid=names(exons_gene_lens),
                                efflen=as.numeric(exons_gene_lens))
    geneid_efflen_gtf1 <- geneid_efflen

2.利用rtracklayer包导入gtf处理

##利用rtracklayer包import导入处理 
    gtf <- as.data.frame(rtracklayer::import("gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz"))
    table(gtf$type)

    exon <- gtf[gtf$type=="exon",
                c("start","end","gene_id")]
    exon_bygeneid <- split(exon,exon$gene_id)   #按照每个geneid所含的exon排序成列表
    
    efflen <- parLapply(cl,exon_bygeneid,function(x){
      tmp <- apply(x,1,function(y){  y[1]:y[2]  }) #输出exon长度值的所有元素            
      length(unique(unlist(tmp))) #去重复并统计exon长度元素的数量
      }) 

    ##转换为dataframe
    geneid_efflen <- data.frame(geneid=names(efflen),
                                efflen=as.numeric(efflen))
    geneid_efflen_gtf2 <- geneid_efflen

所得结果的比较

现在就可以来进行基因有效长度之间的比较啦。
首先看看从gtf文件中获取基因有效长度的两种方法是否有差异。可以发现，仅有极少数efflen有差异，因此这两种方法可以说是几乎没什么差别了：

从gtf文件中获取efflen的比较

再比较一下geneid_efflen_gtf1和geneid_efflen_salmon，发现有一半的efflen是不匹配的？仔细一想这也是可以理解的，因为上游中salmon是对样本中的转录本进行的统计，这说明了样本中有一半的基因并未表达其全部的转录本而已：

salmon和gtf中获取的efflen比较

再将geneid_efflen_gtf1和geneid_efflen_fc进行比较，发现全都能匹配上，这说明featureCounts的Length确实就已经是我们想要的基因有效长度了（即非冗余外显子总长度）：

featureCounts和gtf中获取的efflen比较

总结：

1. 获取基因有效长度的最简便方法是直接从featureCounts或salmon的输出文件中提取。
   但需要注意的是，featureCounts中基因有效长度Length即为基因的非冗余外显子总长度，而salmon中的基因有效长度Length是目标基因的转录本总长度，由于样本中只有部分基因会表达其全部类型的转录本，因此salmon中的转录本总长度会有部分小于非冗余外显子总长度。

2. salmon输出结果中不仅给出了基因有效长度Length（转录本总长度），还给出了经过考虑各种因素矫正后的转录本有效长度EffectiveLength。Salmon官方更推荐使用EffectiveLength进行后续的分析，认为其能更好消除测序时基因长度的影响，它结果中的TPM值也是根据EffectiveLength计算的，后续分析中可以直接采用。

3. 在没有上游原始输出文件的情况下，也可以采取直接从gtf文件中计算的方法，获取每个基因的非冗余外显子总长度得到基因有效长度。还可以保存geneid_efflen便于之后再读取：
   write.csv(geneid_efflen,file = "geneid_efflen_vm25_gencode.csv",row.names = F)

参考资料

基因长度之多少 | 生信菜鸟团 (bio-info-trainee.com)

Counts FPKM RPKM TPM 的转化 - 简书 (jianshu.com)

基因长度并不是end-start - 简书 (jianshu.com)

Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)

Salmon Output File Formats - Salmon 1.8.0 documentation

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年6月场，你的生物信息学第一课